Characterization of the Dynamics of Glass-forming Liquids from the Properties of the 

Potential Energy Landscape 
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We develop a framework for understanding the difference between strong and fragile behavior 
in the dynamics of glass-forming liquids from the properties of the potential energy landscape. 
Our approach is based on a master equation description of the activated jump dynamics among 
the local minima of the potential energy (the so-called inherent structures) that characterize the 
potential energy landscape of the system. We study the dynamics of a small atomic cluster using 
this description as well as molecular dynamics simulations and demonstrate the usefulness of our 
approach for this system. Many of the remarkable features of the complex dynamics of glassy systems 
emerge from the activated dynamics in the potential energy landscape of the atomic cluster. The 
dynamics of the system exhibits typical characteristics of a strong supercooled liquid when the 
system is allowed to explore the full configuration space. This behavior arises because the dynamics 
is dominated by a few lowest-lying minima of the potential energy and the potential energy barriers 
between these minima. When the system is constrained to explore only a limited region of the 
potential energy landscape that excludes the basins of attraction of a few lowest-lying minima, the 
dynamics is found to exhibit the characteristics of a fragile liquid. 



I. INTRODUCTION 

In the supercooled state, glass-forming liquids exhibit 
many fascinating features [1-3] in its dynamic behavior, 
such as multistage, non-exponential decay of fluctuations 
and a rapid growth of relaxation times with decreasing 
temperature. In these aspects, glassy systems challenge 
us with many interesting issues and questions which are 
not well resolved theoretically. 

A popular phenomenological characterization of the 
dynamics of glassy systems, proposed by Angcll, [1-4] 
is the classification of the dynamics as strong or frag- 
ile. It classifies different glass formers on the basis of 
the temperature dependence of their viscosity rj or their 
structural relaxation time r. Quite generally, the rapid 
growth of r with temperature T can be represented by 
a generalized Arrhenius form, t(T) ~ cxp(Ei,(T)/kBT) 
with an effective temperature dependent activation en- 
ergy Eb(T). For strong systems, this activation energy is 
essentially independent of temperature while fragile sys- 
tems exhibit a strong temperature dependence of this 
quantity: it increases as T is decreased. This implies 
that the relaxation mechanism is independent of temper- 
ature for the first type of systems, whereas for the other 
type, it depends on T . Although various measures of the 
extent of fragility exist in literature [5] in terms of the 
slope of the ln(r) vs. 1/T plot, the quantitative distinc- 
tion between strong and fragile liquids at a microscopic 
level is not fully understood yet. 

Following Goldstein [6] , a widely accepted way of look- 
ing at glassy dynamics is to view the dynamical evolution 
of the system in terms of the motion of a state point in 
configuration space, specified by 3 A coordinates for an N 
particles system, over its potential energy surface, often 
referred to as potential energy landscape (PEL) [7]. For 
glass-forming liquids and other disordered systems such 
as spin glasses [8] in general, a generic feature of the PEL 



is the existence of a large number of local minima in it, so 
that the potential energy surface has very rough topog- 
raphy. At low temperature, in the supercooled regime, 
the system visits the neighborhood of a local minimum 
for very long times and makes occasional jumps to other 
minima close to the initial one over the barriers sepa- 
rating them. Following this description, Stillinger and 
Weber [9-11] showed that a useful approach towards the 
understanding of the low-temperature properties of such 
system with a rugged PEL is to divide the configuration 
space into basins of attraction of the local minima and 
then formulate a statistical description in terms of the 
distributions of different properties of these local minima, 
denoted as inherent structures (IS), and their basins of 
attraction. Stillinger also suggested [12] a rationale for 
the difference between strong and fragile liquids based 
on the PEL viewpoint. According to his idea, strong liq- 
uids have a uniformly rough PEL and there is not much 
variation in the values of the energy barriers separating 
different inherent structures. On the contrary, the PEL 
of fragile liquids have non-uniform roughness. At high 
temperatures, the system explores a PEL with nearly 
uniform roughness due to its high kinetic energy, but at 
lower temperatures it explores the deep valleys with very 
different energy barriers, giving rise to a temperature- 
dependent average barrier height. The work of Sastry 
et al. [13] demonstrated the usefulness of this descrip- 
tion, although the understanding of the dynamics in this 
approach still remains qualitative to a large extent. 

After early work [9, 10] on the implementation of a 
description based on inherent structures by supplement- 
ing conventional Molecular Dynamics (MD) simulations 
[14] with regular steepest-descent minimizations of the 
potential energy (called quenches) and thereby sampling 
the inherent structures, a large amount of activity has 
gone into the field in recent years and a variety of meth- 
ods have been suggested to make the survey of the PEL 
more efficient [7]. Once the inherent structures are sam- 
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pled properly, rates of transition between them may also 
be calculated and the dynamics of the system can be 
described by the evolution of the probabilities of occu- 
pying different basins through a master equation [7, 15- 
17]. An important assumption behind such an approach 
is a separation of time scales i.e. the assumption that 
the intra-basin relaxation time Ti ntra and the inter-basin 
relaxation time Ti ntcr are well separated, T; ntra <C Ti n t cr - 
This is indeed a good assumption at low temperatures for 
which the typical barrier height Ef, > kgT. This formula- 
tion does not have the usual limitations arising from the 
requirement of long simulation times, since the master 
equation can be formally solved for all times. However, 
for writing down the rates of transition between two min- 
ima, we need to sample the barriers or transition states 
between them and finding transition states is much more 
demanding computationally than finding the local min- 
ima. However, in recent years, various methods [7] have 
been developed to overcome this difficulty. 

The number of minima, n m i n for model glass formers 
increases exponentially with the number of particles in 
the system [12]. So it is impractical to hope to list all the 
minima even for a system of moderate size, say one con- 
sisting of a few hundred particles. Moreover, the formal 
solution of the master equation requires the diagonaliza- 
tion of matrix of order n min x n min . Hence, the applica- 
bility of this method has been restricted mostly to small 
system sizes till now [7]. Alternatively one can organize 
a set of inherent structures into larger metabasin [19-21] 
and set up a master equation dynamics for transitions be- 
tween the metabasin [21]. Another aim of the metabasin 
construction is to define a space where the dynamics of 
the system point become Markovian, a crucial assump- 
tion behind any master equation based description. The 
Markovian assumption might break down for elementary 
jumps between the basins of inherent structure. However 
assigning the rate of transition between two metabasins 
is not straightforward as one needs to integrate out the 
intra-metabasin dynamics (e.g. jumps between basins in 
the same metabasin) in that case. Also, the construction 
of the metabasins out of inherent structures is generally 
done using somewhat ad-hoc criteria [19-21]. 

Nevertheless, the long-time, low-temperature dynam- 
ics of systems with a small number of particles (say ten to 
hundred) interacting via some model potential has been 
studied very efficiently using this coarse-grained (in time) 
master equation approach since one can make almost ex- 
haustive search of all the minima (a few hundreds to sev- 
eral thousands) and obtain a moderately good number 
of transition states for them [7, 22]. Though small, these 
clusters captures many features of the complex dynam- 
ics [17] observed in larger systems and hence constitute 
a good playground for relating the properties of the PEL 
to the dynamics. 

In the same spirit, we consider here the master equa- 
tion dynamics in a connected network of minima [15], 
where the minima serve as the nodes and transition states 
as the edges of the network [23]. In Section II, we briefly 



review the general formalism [15, 17] for calculating time 
autocorrelation functions of various physical quantities 
and the corresponding relaxation times based on the mas- 
ter equation dynamics in the network of minima. We 
show that a quantitative understanding of fragility can 
be obtained in this framework from an analysis of elemen- 
tary jumps between inherent structures and the effective 
barrier Eb(T) appearing in the temperature dependence 
of the relaxation time can be directly calculated from the 
local properties of the minima and the transition states 
that connect them. We also comment on the breakdown 
of Stokes-Einstein [7] relation observed in many glass for- 
mers from this perspective. 

To test the validity of our results, we study of the equi- 
librium dynamics of a cluster of 13 atoms interacting by 
the Morse potential [24] in Sections III and IV. The PEL 
of this system has been studied in detail in the past [25] 
and a nearly exhaustive list of minima and transition 
states has been obtained. The PEL resembles a funnel, 
in which the minima are organized into pathways of de- 
creasing energy leading to the global minimum. This sys- 
tem has the nice property of having a complex landscape 
that consists of a fairly large but still manageable num- 
ber of minima and also displaying some of the salient 
features of the complicated dynamics of glassy system, 
as we report in Section IV. By restricting the system to 
sample certain parts of the PEL excluding a few lowest- 
lying minima, we are able to show that the dynamics 
of the system in this restricted part of the configuration 
space exhibits characteristic behavior of fragile systems, 
i.e. Ef,(T) is perceptibly temperature dependent, whereas 
the dynamics in the full PEL exhibits strong features. 
We have also carried out MD simulation for the 13-atom 
Morse cluster. Using simulations in which the MD tra- 
jectories are confined [26, 27] in appropriately restricted 
parts of the PEL, we are able to substantiate the conclu- 
sions of the network model calculation. The results from 
MD simulations are discussed and compared with those 
obtained from the network model in Section V. Some of 
the technical details of the network model calculations 
and restricted MD simulations are described in the Ap- 
pendices A, B and C. 

II. MASTER EQUATION FOR JUMP 
DYNAMICS BETWEEN INHERENT 
STRUCTURES 

The first step in exploring the landscape is to find 
the configuration corresponding to the local minima and 
the transition states. These are the stationary points or 
saddles of the potential energy function V(ri, r 2 , rjv) 5 
characterized by W = and number of negative eigen- 
values of the Hessian matrix H = V?/ = ® Y h , rf 

being the a-th coordinate of the i-th particle. The eigen- 
vector corresponding to each negative eigenvalue of the 
Hessian matrix at a stationary point signifies an unsta- 
ble direction and stationary points can be indexed by 
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the number of such unstable directions. For instance, at 
a local minimum, there is no negative eigenvalue and it 
can be denoted as saddle of index 0, a transition state as 
saddle of index 1 or first-order saddle and similarly there 
can be higher order saddle points having index running 
from 2 to the dimension of the configuration space. One 
can use steepest descent minimization [28] for finding the 
minima and transition states can be located efficiently 
using the eigenvector following method [29, 30]. Once 
the minima and transition states that connect them are 
known, we can arrange them by designating as nodes and 
edges, respectively, of a network. We describe this pro- 
cedure in greater detail later for the specific case of a 
13-atom Morse cluster (see also Appendix B). 

In the model of a connected network of potential en- 
ergy minima, the master equation for the jump dynamics 
can be written as 

^ = $> oc P c (*;Mo). (1) 

c 

P a (t;b,to) is the probability that the system is at min- 
imum a at time t, if it was at a minimum b at 
time to and a runs from f to n m i n , the total num- 
ber of minima in the network. The off-diagonal ele- 
ments of the matrix W are the transition rates and, as 
usual, the diagonal elements are fixed by the condition 
Ea P a(i;Mo) = 1 implying J2 a W ac = 0. In order to 
obtain an asymptotic behavior (in the long time limit) 
that agrees with the Boltzmann distribution, the occupa- 
tion probability must satisfy lim t _ K3C) P a (i; b; to) = P® = 
Z- 1 (Det(H a ))- 1 / 2 exp(-K/T). Here, Z is such that 
J2 a P^ = 1 an d the pre-exponential factor follows from 
a harmonic approximation for the partition function in 
the basin of each minimum (we take the Boltzmann con- 
stant kg = 1). The matrix H a is the Hessian matrix 
for the o-th minimum. As W ought to satisfy the de- 
tailed balance relation W a bP° = Wb a Pa, in numerical 
calculation it is more convenient to express the solution 
in terms of the eigenvectors of the real symmetric ma- 
trix W a b = W a b(P b °/P°) 1 / 2 . Finally, one can formally 
obtain the solution of the master Equation (f) for all 

time, i.e. P a (t;b,t ) = (i*/i?)^E«e£ n) 4 n) e A » ( '-*°>. 
Here e^™) are the eigenvectros of W corresponding to the 
eigenvalue X n (1 < n < n min ). The matrix W has one 
zero eigenvalue corresponding to the equilibrium distri- 
bution and all other eigenvalues are negative. We shall 
follow the convention of arranging A n 's in descending or- 
der (ascending order in their absolute values) starting 
from Ai = 0. 

The model is well-defined once we give an appropriate 
expression for the transition or hopping matrix W be- 
tween the nodes of the network of minima, Treating the 
problem as a Markovian Brownian multi-dimensional mo- 
tion in the over-damped friction regime [31-33], we can 
write the transition rates between a directly connected 
or nearest-neighbor pair of minima, < ab >, from b to a 



over the saddle s as 

_ fl. a , nt / Det(H a ) ^ ( ) 

Here, Cj s . a b is the down frequency at the saddle point 
i.e. a) 2 ab — A* b , A* b being the magnitude of negative 
eigenvalue of the Hessian matrix H* b at the transition 
state, n is the friction constant that sets the time scale 
(we take jj, = 1 in all our calculation henceforth), and 
Vb, V Q s b are, respectively, the potential energies at the 
minimum b and the saddle point s between minima a and 
b. If there are multiple barriers connecting b and a then 
the total transition rate from b to a, W a b is obtained by 
summing over all the barriers i.e. W a b = J2se<ab> ^ab- 

A. Correlation function 

To study the dynamics and calculate relevant relax- 
ation times in the network model, one needs to de- 
fine the equilibrium time autocorrelation function of 
some physical quantity, say c/>(r(t)), a generic observ- 
able which depends on the collective coordinate r(t) = 
{i"i(i), r 2 (i), i"at(£)} at times t. In this language, the 
time autocorrelation function can be written [15, 16] as 

(c/>(r(t))cl>(r(to))) = m,to)) = Ct(t,t ) 

b a 

= C° + £ e*"^ E *a 6 (P o °P h ) 1/2 4 n) 4 n) - (3) 

n>2 a, b 

Here $ a6 = $(r a ,r 6 ) = 0(r a )0(r 6 ) and C° = 

Ea^abW- AISO, CV(f = t ,t ) = Ea^aaPS and 

lirrif-j.oo C^(t,to) = are short-time and long-time lim- 
its of the correlation function C^,(t,to), respectively. 

Once the correlation function is calculated, the relax- 
ation time can be estimated by assuming a pure Dcbyc 
(single exponential) relaxation, such that C$(t — to) = 
C<t,(t,t ) ~ C° = C^(t ) exp (-(t - to)/r|) and evaluat- 
ing the area under the resulting curve from Eq.(3), i.e. 

1 I pO p(m/2rh An) (n) 

r e _ Ln>2 IT l^a,b\ a b ) ^ab^a e b ,^ 

L-m>2,a,b\ a b I v ab&a C b 

This way of defining the relaxation time relies on the as- 
sumption that the decay of C$(t — 1 ) is well described by 
a single exponential, but in glassy systems the temporal 
decay of correlation often follows a profile that is more 
complex than a simple exponential. The precise form of 
the decay is usually not known, although it can be fitted 
with the empirical Kohlrausch- Williams- Watts (KWW) 
or stretched exponential function [34, 35] 

c4t) = cy^ (5) 
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in many cases. Here r$ is a measure of the relaxation 
time and /3 is the KWW (or stretching) exponent. The 
above form [Eq.(5)] and other typical features of glassy 
dynamics such as the decoupling of transport coeffi- 
cients (i.e. the violation of the Stokes-Einstein relation 
rj (x T/D, r) and D being the viscosity and diffusion coef- 
ficient, respectively) can be rationalized from the hypoth- 
esis of the existence of a distribution of time scales that is 
not sharply peaked at a particular value. A distribution 
of time scales may arise e.g. from dynamical heterogeneity 
that describes spatial variations of the local relaxational 
dynamics (i.e. the fact that different parts of the sample 
may have different relaxation times). Given a distribu- 
tion of relaxation times, p^(t), the temporal decay of a 
typical time autocorrelation function, reflecting the ef- 
fects of all these relaxation processes, is given by 

/>oo 

C^t) = / e-Wp*(T)dr, (6) 
Jo 

and one can easily show that this results in a stretched 
exponential form for an appropriate distribution p^(t). 
In the present master equation framework, a distribution 
of relaxation time is naturally provided by the different 
modes n = 2, ...,n m ; n and one can identify the contri- 
bution or weight [similar to p (t) in Eq.(6)] of the n-th 
mode as 

(n) _ 2^a,b\ r a r b! V ab e a e b 

W<i> ~ V fpopou/2$ , f W P W ! 

2^n>2,a,iA a b > ^ab^a e fe 

where r„ = |A„| _1 is the relaxation time corresponding 
to the n-th mode. A distribution of relaxation times 
in the jump dynamics on the PEL appears due to a het- 
erogeneous distribution of barrier heights contributing to 
different relaxation channel n. Though the origin of non- 
exponential relaxation due to a broad distribution of re- 
laxation times in the master equation based framework 
looks superficially similar to that due to dynamical het- 
erogeneity mentioned above, the precise correspondence 
between the heterogeneity of barrier heights in the con- 
figuration space and dynamical heterogeneity in the real 
space is not very clear. 

As shown in Appendix A, one can utilize Eq.(7) to 
directly calculate the temperature dependent activation 
energy Eb(T) [or, more precisely E b (T)] to a good ap- 
proximation and a quantitative measure, namely W^ b , of 
the participation of individual barriers (or the pairs of 
minima that are connected by the barrier) in the relax- 
ation process. We have tested these results [Eqs. (A2a) 
and (A2b) in Appendix A] for the case of 13-atom Morse 
cluster in Section IV. 



B. Structural relaxation 

Since we are mainly interested in structural relaxation 
of the system as its state point explores different regions 



of the PEL, we look for the time autocorrelation function 
of quantities related to the configurations of the min- 
ima. One such quantity whose autocorrelation function 
and related relaxation time is of much interest in glass 
physics is the off-diagonal microscopic stress tensor [36], 
i.e. = E^vf - Y Ji<j V'{r ij ){rt :j r ij /n j ) {a ^ /3). 
Here, vf and are the a-th components of the velocity 
of the z-th particle and r.y = r i — rj , respectively, while 

V'{r) = dV Qr ■ Since the velocity v.; is not defined in the 
model, following reference [15], we neglect the kinetic en- 
ergy term and work with the following quantity 

= -J^V'in/-^. (8) 

Consequently, we can calculate the stress-stress auto- 
correlation function C a {t) = (1/3) E a </3( cr " /3 W^W) 
from Eq.(3) where we need to insert $ a {, — 

J2a<p <J a^ (T b ■ The shear viscosity, rj oc T _1 dtC a (t), 
is related to the stress autocorrelation function. We 
compute another correlation function, termed as over- 
lap function C$(t), by setting $ a (, = $(r(t) <E £> a ,r(0) <E 
13b) = S a b in Eq.(3) to meaningfully compare the pre- 
diction of the network model with MD (see Section V). 
The basin of attraction B a of the a-th minimum is the 
set of state points that flow to the a-th minimum un- 
der steepest descent minimization and 8 a b is the usual 
Kroneckcr delta function. The correlation function Cg(t) 
decays solely due to transitions between inherent struc- 
tures and it can be calculated from MD simulation as we 
show in Section V. 



C. Diffusion constant and waiting times 

The diffusion constant D can be calculated from 
the mean square displacement (R 2 (t)) using R 2 b = 

N^Etil^ " rf| 2 in place of $ ab in Eq.(3) 
(r^ denotes the position of the i-th particle in the 
a-th minimum) and invoking the definition D = 

linw^^ = lim^i^. This def- 
inition would yield D = here as D = 

^^ooEaAnUPSPS^Kb^e^e^nX^t) = 0. 

In other words, (R 2 {t)) ~ H*) and lim^oo a(t) = 0. 
Instead, we can define the diffusive regime as the time 
window over which a(t) ~ 1 and then extract D from 
the slope of the (R 2 (t)) vs. t plot in this region. This 
may not be useful for confined systems like the atomic 
cluster considered here, since such a diffusive regime may 
be very short or absent altogether, rendering D to be an 
ill-defined quantity. Nevertheless, in this model we can 
define D from the initial slope of the (R 2 (t)) vs. t curve 
as the usual ballistic regime ((R 2 (t)) oc t 2 ), seen for in- 
stance at very short times in MD simulations, is absent 
by definition. So, we may define the diffusion constant 
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FIG. 1: The essential PEL details, that go into the calculation of transition rates [Eq.(2)], are shown. Panel a: Histogram of 
IS energies. Panel b: Dependence of the overall curvature of a minimum on its energy. Panel c: Heights of the lowest energy 
barriers for going from minimum a (with potential energy V a ) to minimum b (with potential energy V&). 



a.s 



D 



lim 

i-s-0 



d(R 2 (t)) 
dt 



= Y,(PaP° b ) 1/2 Rl b W ab . (9) 
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The waiting time in the basin of a minimum is the 
amount of time the system spends between an entry into 
and the subsequent exit from the basin (i.e. during a sin- 
gle visit to the basin). The average of this quantity for 
the a-th minimum, (T w (V a )), can be calculated both from 
MD simulations and in the network model. In the mas- 
ter equation based model it is straightforward to write 
(r w (V a )) as 

(r w (V a )) ee r a = -^—. (10) 

V* aa 

On the other hand, if we consider an crgodic MD trajec- 
tory, then the amount of time the system spends in B a is 
t a oc P°. Hence the total number of visits to B a over the 
full trajectory can be written as v a ee (t a /T a ) oc (Pa/ T a)- 
Finally, the mean waiting time averaged over the whole 
landscape would be 

Ea T « U a _ 1 



The above quantity can be related to D [Eq.(9)] 
if R^ b ~ R 2 for all pairs of minima, as D ~ 

R 2 EaM P a P b y /2W «>> = -^EaPSWaa = R 2 /t w . 

Hence for this special case the hopping between the 
basins becomes a random walk with a distribution of 
waiting times [37]. 

III. CONSTRUCTION OF THE NETWORK FOR 
A 13-ATOM MORSE CLUSTER 

The Morse potential [24] can be written in the form 



g p(l-ry/r e ) r e p(l-rij/r e ) 



2]e,(12) 



where is the distance between atoms i and j, e and r e 
are the dimer well depth and the equilibrium bond length, 
respectively - they simply scale the PEL without affecting 
its topology and can conveniently be set to unity and used 
as the units of energy and distance. The parameter p is a 
dimensionless quantity that determines the range of the 
potential, with low values corresponding to long range. 
We have taken p = 4. This potential is widely used to 
model inter-atomic interactions in small atomic clusters 
or molecules [7] . The reduced unit of time can be set to 
(mfg/e) 1 / 2 , m being the atomic mass. 

We follow more-or-less the same procedure (described 
briefly in Appendix B) as that mentioned in Ref.[25] for 
building the network of minima and transition states. 
The network that we obtain consists of 138 minima and 
230 transition states connecting them. We have not en- 
forced any confining potential to prevent the cluster from 
melting - in the temperature range of our interest, the 
particles are confined due to interactions among them- 
selves (we have discarded minima with maximum inter- 
particle separation more than 2.5). The minima and the 
transition states around a particular minimum are iden- 
tified by the values of their potential energy 

In Fig.l a, the distribution of inherent structure en- 
ergies is shown. We index the inherent structures in as- 
cending order of energy. The lowest lying minimum (the 
global minimum) is at V a = —46.635 and after a substan- 
tial gap there are three minima at V a — —43.5; after that, 
another perceptible gap is present and the rest of the 
minima are closely spaced for V a ~ — 42.5. Henceforth, 
we denote by Aff the full network and by J\f r the net- 
work with the four lowest-lying minima removed. When 
we remove a particular minimum, all its edges and min- 
ima that are connected to the rest of the network solely 
through this particular minimum get deducted from the 
network as well. As a result, Af r contains 52 minima and 
70 transition states. 

The overall curvature at a minimum, C is obtained 
from the determinant of the Hessian matrix at the min- 
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imum i.e. C = Dct(H). Here the deeper basins are nar- 
rower, as is evident from Fig.l b. Also the number of 
barriers and barrier heights increase with decreasing IS 
energies (Fig.l c). 

In the next section we describe the results of the master 
equation based calculations carried out with the networks 
Aff and M r - While the relaxation dynamics in Mf is 
quite similar to that of strong glass formers, the dynamics 
in M r shows strong resemblance to that of fragile ones. 
Although Mr is a part of A//, the relaxation of the system 
restricted in this part of the configuration space becomes 
qualitatively very different from the global dynamics in 

M f . 

The procedure of searching for transition states (Ap- 
pendix B) starting at a minimum and moving along di- 
rections of successively larger eigenvalues of the Hessian 
matrix is not very efficient [22] and most of the time, 
one ends up at the same transition state. Hence we get 
only 3 barriers around a minimum on the average. Al- 
though the details of the connectivity of the network may 
matter for finer details, our main objective in this work 
has been to realize the typical characteristics of glassy 
dynamics with a minimal set of minima and transition 
states. Hence we have not paid much attention to build- 
ing the network very accurately so as to represent the 
actual system. Nevertheless, the network model seems 
to capture the main features of the complex long-time 
dynamics quite accurately, as verified through MD sim- 
ulations (Section V). 



IV. STRONG AND FRAGILE BEHAVIOR IN 
THE NETWORK MODEL 

Here we describe the results for the stress-stress auto- 
correlation function C a (t), as described in Sec.IIB, for 
the Morse cluster. 

Fig. 2 a shows C a (t) for the network Mf at four different 
temperatures. The decay of the correlation function is 
very well described by a single exponential [the stretching 
exponent of Eq.(5), /3 ~ 1] over the entire temperature 
range, as is evident from the inset of Fig. 2 c. The origin of 
this simple Debye-like relaxation can be attributed to the 
presence of a single relaxation mode with a large weight 
w„ [Eq.(7)] as shown in Fig. 2 b. The corresponding 
relaxation time r CT extracted by fitting C a (t) with Eq.(5) 
follows an Arrhenius temperature dependence (Fig. 2 c), 
i.e. 



Tcr ex P ( Y 



(13) 



In Mf, the dynamics is mainly governed by the barri- 
ers connecting the global minima at V a ~ —46.6 with 
the next three lowest lying minima at V a ~ —43.5 (see 
Section III). 

The Arrhenius fit to r a vs. T curve yields an effective 
activation barrier E£ = 0.24 which one can easily iden- 
tify with the barrier for going from one of the minima 
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FIG. 2: Stress autocorrelation function for the network Mj 
(Section III). Panel a: C CT (t) [normalized by dividing with 
Cct(0)] is shown for various T. Panel b: Decay of C a (t) is 
determined by one dominant relaxation mode n over the en- 
tire temperature range of interest, as is evident from the plot 
of io?° vs. |A»| -1 [Eq.(7)] for T = 0.10, 0.25. Panel c: Ar- 
rhenius plot for T a (T) [extracted by fitting C a (t) with the 
KWW form of Eq.(5)j i.e. lnr CT vs. 1/T. The effective bar- 
rier E% is obtained by fitting the data to the Arrhenius form 
[Eq.(13)]. The stretching exponent /3 ~ 1 (inset) confirms the 
simple exponential nature of the decay. Panel d: The esti- 
mate of Eq.(A2a) (Appendix A) for E% agrees very well with 
the value of E% = 0.24 obtained from the Arrhenius fit shown 
in panel c. 



at V a — —43.5 to the global minimum. The estimate 
for El (Fig.2 d) obtained from Eq.(A2a) (Appendix A) 
agrees very well with the above value. In contrast, we 
find that the barrier for going from the global minimum 
to the next lowest lying minimum determines the effec- 
tive activation barrier that appears in the T-dependcncc 
of the mean waiting time t w [Eq.(ll)]. This is a natural 
consequence of the fact that the global minimum, being 
much lower in potential energy with respect to the other 
minima in this case, possesses almost all the Boltzmann 
weight. As a result, the relaxation to equilibrium (decay 
of correlation) is entirely dictated by the relaxation paths 
from other parts of the PEL to the global minimum. The 
mean waiting time t w , on the other hand, is decided by 
the escapes from the global minimum over the barriers 
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surrounding it as the system spends most of the time in 
the basin of the global minimum. 

The above observations suggest a trivial route for re- 
alizing strong behavior, namely dynamics governed by 
a fixed set of barriers surrounding a very deep inherent 
structure (or a set of inherent structures with very simi- 
lar potential energies in a more general case) possessing 
most of the Boltzmann occupation probability. Keeping 
this fact in mind we construct the network M r by remov- 
ing a few deep minima so that a larger number of minima 
figure in the relaxation to equilibrium due to compara- 
ble Boltzmann weights in the activated regime and many 
different barriers contribute to the relaxation process. 

Fig. 3 a exhibits C a (t) calculated for Af r . We observe a 
two-stage, non-exponential relaxation in this case. This, 

again, can be understood from the values of (Fig. 3 
b). The profile of C a (t) is well fitted with a sum of two 

stretched exponentials, i.e. C a (t) = C^\t) +C„\t) = 

We plot the 



C«(0> 



-(*/• 



> ft + Ci 2) (0)e-(*M 2) )' 32 



( 2 ) 

longer of the two relaxation times, t& , in Fig. 3 c along 
with the exponent /?2 in the inset. The deviation from 
simple exponential behavior is evident. The temperature 

(2) 

dependence of t£ exhibits marked deviation from the 
simple Arrhenius behavior. Rather, the Vogel-Fulcher- 
Tammann (VFT) form [38-40], 



r°exp 



T - Tk 



(14) 



frequently used for fragile glass formers, yields a much 

(2) 

better representation of the data for ts- . The fragile 
nature of the hnv vs. 1/T curve is well reproduced in 
Fig. 3 d by the effective temperature dependent barrier 
E°{T) calculated from Eq.(A2a). 

A useful visualization and understanding of the ob- 
served strong behavior for the PEL A/"/ and fragile be- 
havior for M r can be achieved by looking at the quantity 
W a b defined in Eq.(A2b) and the subsequent paragraph 
of Appendix A. This quantity, W a b, can be thought of 
as the weight with which the edges between the nodes a 
and &, or in other words, the elementary jumps over the 
barriers connecting minima a and b appears in the re- 
laxation process through various relaxation channels or 
modes n. In Figs. 4 a, b and c we have shown the color 
maps of W a b for Mf at three temperatures. In these 
plots, the coordinates (V a ,Vb) represent the barrier be- 
tween minima a and b and the color corresponds to the 
value of W a b at (V a , 14) (we have used a small broaden- 
ing for the purpose of visualization) . It is clear from the 
plots that only a few barriers surrounding the global min- 
imum contribute substantially to the relaxation process 
and over the entire temperature range, the peaks remain 
at nearly the same positions. On the contrary for M r , in 
Figs. 4 d, e and f, many barriers figure in the activated 
relaxation and the picture changes substantially with in- 
creasing temperature, as more and more barriers start to 
play a role in the relaxation process. This feature can be 
considered as the trademark of fragile dynamics. 
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FIG. 3: Stress autocorrelation function in Mr- Panel a: C CT (t) 
is shown for five T values. Panel b: Fit to a sum of two 
stretched exponentials (Eq.(5), see the text for details) at 
T — 0.08. The quantity exhibits the presence of (well- 
separated) multiple timescales even at the very low tempera- 
ture T = 0.08. Panel c: The deviation from the simple Arrhe- 
nius form is evident from the plot of hiri- 2 ^ vs. 1/T. Fits to 
both the Arrhenius form [Eq.(13)] and the VFT form [Eq.(14)] 
are shown. It is clear that the VFT form provides a good fit. 
Inset: The exponent /3a is much less than 1, specially at low 
temperatures, showing the non-exponential nature of the re- 
laxation. Panel d, The effective barrier E%(T) obtained from 
Eq.(A2a) agrees reasonably well with the estimate deduced 
from the Arrhenius and VFT fits. For instance, the VFT fit 
to exp(E£/T) vs. 1/T yields values for the parameters B CT 
and Tq [Eq.(14)] that are similar to those obtained in panel 
c. 



We have observed similar characteristics of strong and 
fragile dynamics for the relaxation time t$ associated 
with the overlap function Cs(t), defined in Section II B, 
obtained from both the network model calculation and 
MD simulations. We report these results in the next sec- 
tion. 
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FIG. 4: Color maps (see the color bars) for W a b [Eq.(A2b)] in the (V a , Vb) plane for A// at three temperatures, T — 0.15 (panel 
a), T = 0.20 (panel b), and T = 0.25 (panel c), and for Af r at T = 0.15 (panel d), T = 0.22 (panel e), and T = 0.28 (panel f). 



V. COMPARISON OF THE RESULTS OF THE 
NETWORK MODEL WITH THOSE OF MD 

We have carried out MD simulations [14] for the 
cluster of 13 particles interacting via the Morse poten- 
tial [Eq.(12)]. Initial long MD runs along with regular 
quenching provide the starting data base of minima and 
then we follow the procedure of Appendix B to improve 
the list as well as to construct the network of minima 
and transition states. The main motivation of our MD 
study is to check to what extent the dynamics of the 
systems is captured by the network model by compar- 
ing the results of the network model calculation with 
those for the real dynamics (i.e. Newton's equation of 
motion). We find that the main results of the network 
model calculation are supported both qualitatively and 
quantitatively by the MD results. As already discussed 
in the introduction, the model assumes that there is a 
clear separation between the time scales of local vibra- 
tions and activated jumps. This assumption is bound to 
be valid at very low temperature where the barriers are 
much higher than the temperature. However, the barrier 
heights have a wide distribution (Fig.l c) and hence the 
assumption of well separated time scales may be invali- 
dated at different temperatures for different parts of the 
landscape. Presumably, as the temperature is increased, 
the dynamics crosses over from a landscape dominated 
regime at low temperatures to a high-temperature regime 



where the system does not see the details of the PEL due 
to its large thermal energy and the dynamics is no longer 
dominated by activation processes. 

We have calculated the stress autocorrelation function 
C a (t) and the overlap function C$(t) (Sec. II B) from MD 
simulations at different temperatures. Also, we have esti- 
mated the mean waiting times and waiting time distribu- 
tions. We find that due to intra-basin relaxation present 
in MD, the stress autocorrelation function C a (t) decays 
within the vibrational time periods of the local minima 
much before the inter-basin jump dynamics comes into 
play. As a result, C a (t) is not useful for comparison 
with the network model results. Hence we concentrate on 
Cs(t) and the related relaxation time since by definition 
the decay of C$(t) is entirely determined by inter-basin 
transitions. For calculating this quantity as well as the 
waiting times for different minima, we construct a discon- 
tinuous trajectory in terms of configurations of inherent 
structures from the MD trajectory and track the transi- 
tions along it by the interval bisection method described 
in reference [19] [see Appendix C]. 

For comparing the results obtained for Af r (Sec. Ill), 
additionally, we need to confine the MD trajectories in 
a restricted part of the configuration space so that the 
basins of the four lowest-lying minima are not visited 
(see Appendix C). This goal is achieved by applying a 
procedure similar to that described in Refs [26, 27]. The 
application of this method in conjunction with the in- 
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FIG. 5: Overlap function Cg(t) for the network A//. Panel a: 
(C 6 (t) - Cs)/(C s (0) - Cg) obtained from the network model 
for four temperatures. Panel b: MD results for Cs(t). Panel 
c: Arrhenius plot for ts(T) deduced in a from KWW fits 
[Eq.(5)]. The top and bottom insets show, respectively, the 
estimates of E^(T) from Eq.(A2a) and the stretching expo- 
nent /3 for the KWW fit to Cs(t) in a. Panel d: Estimates 
of t$ can be obtained accurately only for a narrow tempera- 
ture range in MD, as discussed in the text. The exponent ft, 
extracted from b, is shown in the inset. 



terval bisection method makes the MD runs very time 
consuming and hence it is difficult to get very long MD 
trajectories (typically 10 8 or 10 9 MD steps with the MD 
step length St = 0.001). So we have averaged over many 
parallel runs with smaller number of MD steps (2 x 10 7 for 
A// and 2 x 10 6 to 2 x 10 7 for Af r , depending on the tem- 
perature) starting from different parts of the landscape 
i.e. different initial conditions (different IS configurations 
are taken from the already existing data base) . 

We compare the results for Cs (t) obtained through the 
network model for A/"/ and Af r with those obtained via 
isokinetic MD [14] in Figs. 5 and 6, respectively. Fig. 5 
a shows the network model results for Cs in A//. Cor- 
responding MD results are shown in Fig. 5 b. The Ar- 
rhenius plots for associated relaxation times are shown 
in Figs. 5 c and d. The strong behavior is evident and 
decays of Cs{t) in both the cases are close to exponen- 
tial ((3 ~ 0.9) as exhibited in the insets of these figure 
panels. A good estimate of the effective activation en- 



FIG. 6: Overlap function Cs(t) for the network M r - Panel 
a: Cs(t) calculated in the network model exhibits multi-step 
decay of correlation. Panel b: Cg(t) as obtained by confining 
the MD trajectory in a part of the configuration space that 
excludes the basins of attraction of the first four lowest-lying 
inherent structures. Panel c: The intermediate relaxation 
time rj 2 ' obtained by fitting the data for Cs (t) in a with the 

sum of three stretched exponentials. Here the mr] 2 ' vs. 1/T 
plot shows a small bending indicating deviations from the 
Arrhenius T-dependence. The VFT fit is also shown, the up- 
per and lower insets show, respectively, the E£(T) calculated 
from Eq.(A2a) and the exponent /?2 from the fits to Cs(t) at 
different temperatures, d, The results for t$(T) obtained in 
MD simulations agree with the fragile behavior observed in 
the network model. The semiquantitative agreement can be 
verified by comparing the VFT fits for ts(T) in the network 
model and in MD. 



ergy is once again obtained from Eq.(A2a) (inset of Fig. 5 
c). Here the effective barrier ~ 0.25 is comparable 
to T in the temperature range of interest and hence the 
associated time scale for relaxation is quite short. Even 
then, an accurate calculation of the correlation function 
C$(t) and the relaxation time ts in MD is quite diffi- 
cult because the system gets trapped in the deep global 
minimum most of the time and a proper sampling of the 
relevant relaxation paths connecting the higher minima 
to the global one becomes hugely time consuming and 
increasingly difficult with decreasing temperature. As a 
result we can estimate ts only for a limited range of tern- 
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peratures, as shown in Fig.5d. Also for temperatures 
higher than T ~ E b ~ 0.25 the network model deal- 
ing only with the activated processes becomes unreliable, 
rendering a comparison with MD results inappropriate. 

For the network Af r , Cs(t) exhibits a three-stage re- 
laxation profile (Fig. 6 a) owing to the presence of well- 
separated time scales similar to the case of C a (t) (Fig. 3 
b). Such a clear cut separation of timescale is not ob- 
served in the corresponding microcanonical MD [14] re- 
sults for Cs(t) shown in Fig. 6 b. although a faint sig- 
nature of multi-step relaxation can be seen at low tem- 
peratures. The decays of Cs(t) at various temperatures, 
shown in Fig. 6 a, are fitted to a sum of three stretched 
exponentials [Eq.(5)]. The temperature dependence of 

(2) 

the intermediate relaxation time rl (T) shows deviations 
from the Arrhenius form [Eq.(13)] in Fig. 6 c, although 
much less in extent than that for (T) (Fig. 3 c). The 
latter fact is indicative of the importance of the quan- 
tity <j> & of Eq.(3), i.e. the quantity whose autocorrela- 
tion function is used to calculate the relaxation time, in 
determining the degree of fragility obtained from the tem- 
perature dependence of the relaxation time. As discussed 
above, the distribution and temperature dependence of 
the quantities W% b , which depend on the quantity cf>(r) 
whose autocorrelation function is used to calculate the re- 
laxation time (see Appendix A), play an important role 
in determining the degree of fragility. We have checked 
that plots of W a b for the overlap function, similar to those 
shown in Fig. 4, exhibit a less pronounced dependence on 
the temperature. This is consistent with the observation 
that the temperature dependence of the relaxation time 
extracted from the decay of the overlap function exhibits 
smaller deviations from the Arrhenius form (lower degree 
of fragility) in comparison to that for the relaxation time 
obtained from the stress autocorrelation function. The 
Arrhenius fit parameter E b and VFT fit parameters B$ 
and T a are shown in Fig. 6 c. These agrees reasonably 
well with the results of similar fits done for t$ estimated 
through MD simulations (Fig. 6 d). 



VI. CONCLUSION 

Our work shows that the master equation based ap- 
proach, first proposed by Angelani et al. [15], for the 
landscape dominated activated dynamics on a network 
of minima and transition states is capable of addressing 
many important issues and challenges related to glassy 
dynamics. Our results for the full network A/ are simi- 
lar to those of Rcfs.[15, 16] where strong behavior in the 
dynamics, arising from the dominance of the global mini- 
mum and the barriers surrounding it, was observed. Our 
study of the restricted network M r leads to the impor- 
tant result that the master equation approach may lead 
to fragile dynamic behavior if many inherent structures 
and the barriers between pairs of them are involved in 
the relaxation process. This conclusion is confirmed by 



our MD simulations. Our study, thus, provides valuable 
insights into the origin of fragile dynamic behavior in 
glass-forming liquids. It also illustrates the usefulness of 
the master equation approach in studies of glassy dynam- 
ics. Due to computational constraints, the applicability 
of this approach has been limited to small system sizes 
till now. A related approach [41] that differs from the 
one considered here in the details of its implementation 
seems to be a promising one for studying larger systems. 
Another interesting future direction for studying moder- 
ately large systems might be to extend this framework for 
the dynamics in the metabasin space to obtain quantities 
relevant for the characterization of glassy dynamics. 
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Appendix A: Effective activation energy 

As mentioned in Section II A, t„ = |A„| _1 is the re- 
laxation time corresponding to the n-th mode (n > 2) 



and w)[ lJ [Eq.(7)] provides the distribution of such re- 



in) 

laxation times. We assume that the relaxation time for 
each individual mode is governed by an effective barrier 



(n) 



i.e. t„ ex exp (E b 



(n) 



/T). Hence we can estimate E 



00 



approximately from the local slope of the lnr„ vs. 1/T 
curve as follows 



E, 
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To obtain the expression (Ala) for E b n ^ we have rewritten 

A„ as J2 a b Wabea 1 * 1 and while taking the derivative of 
A„ with respect to T neglected the temperature depen- 
dence of the eigenvector e(™) assuming it to be weakly 
temperature dependent. Wc find this to be valid for the 
case of the 13-atom Morse cluster and the above approx- 
imation for E b seems to work, as we have reported in 
Sections IV and V. The overall effective activation energy 
E b (T) is obtained by summing over the contributions of 
all the modes appearing with the weights w^' [Eq.(7)], 



(Ala) 
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En>2 E b 

Et 



(n) (n) 



This can be recast as 



= E 



b> nj <ab> 



(A2a) 



<ab> 
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here, 



yy <ab> 



E 



<ab> 



/ , * v s,<ab> 
s 

Z-is c 's,<a6> 
t—ts ' v s,<a 



(A2b) 
(A2c) 



with 



s,<a6> 



£s,< ab > = e iAnrM n) te^irJ -z4?) 2 



n>2 



The effective barrier E%(T) obtained from Eq.(A2a) com- 
pares quite well with the barrier extracted from the In- 
dependence of the relaxation time estimated by fitting 
the KWW form [Eq.(5)] to the correlation function C^,(t), 
computed using Eq.(3) for the Morse cluster. The quan- 



tity VC 



W< ab> /2 can be a good measure of 



the relative importance of a pair of minima or an elemen- 
tary jump in the relaxation process as demonstrated in 
Section IV. 



Appendix B: Construction of the network of minima 
and transition states 



As mentioned in Sec. Ill, we follow the procedure simi- 
lar to that in Rcf. [25] for building the network model for 
the 13-atom Morse cluster. We have used the OPTIM 
package [30], developed by D. J. Wales and co-workers, 
to search for minima and transition states as follows 

1. Do long MD simulations and steepest descent 
quenches in regular intervals to reach nearby min- 
ima along the MD trajectories. Distinction between 
configurations of minima is done in terms of their 
potential energy values. This provides us with an 
initial data base or list of minima. 

2. Starting from each of these minima, one by one, 

(a) search for a transition state along the eigen- 
vector with the lowest eigenvalue. 

(b) After reaching a transition state, do steepest 
descent minimization starting parallel to the 
transition vector, i.e. eigenvector correspond- 
ing to the negative eigenvalue, at the transi- 
tion state to arrive at the minima connected 
directly to initial one (sometimes, following 
the eigenvector from a minimum we obtain 
states which are not connected to it by steep- 
est descent - we discard these states). This 
establishes one edge, constituted of a pair of 
minima connected by a transition state, of the 
network. 



(c) Repeat (a) beginning anti-parallel to the 
eigenvector with the lowest eigenvalue and 
then successively in both directions along 
eigenvectors with ascending eigenvalues until 
all the eigendirections are considered at the 
starting minimum. 

(d) If some new minima, not in the starting list, 
are found in this process we append them to 
the data base and modify the existing list of 



3. Repeat 2 until all the minima in the data base are 
searched. 



Appendix C: Interval bisection and confinement of 
MD trajectories in a specified part of the PEL 

Starting from a MD trajectory, r(t) = 
(ri(£),r 2 (£), ...,rjv(t)), we construct a trajectory 
r°(t) = (r?(i),r!](i),...,r^(i)) in terms of the inherent 
structure configurations by doing steepest descent 
minimization at regular intervals [19]. The straightfor- 
ward way, though computationally impractical, is to 
do minimization at every time step. Also there may 
be many back and forth jumps between neighboring 
minima giving rise to events that do not affect the 
long time relaxation which we are mainly interested in. 
Bearing this fact in mind, we do quenches for equidistant 
points, t n = nAt (n being an integer) to get r°(t n ) with 
At = 100 - 200 MD steps (around 1/10-th of the typical 
vibrational period at a minimum for the Morse cluster). 
During a MD run, starting from initial point (n — 0), 

1. Save the trajectory r(t) between = t n and tf = 
t n +i and then follow the steps described below pro- 
vided r°(tj) 7^ v °(tf) (actually we compare V(r°) 
instead of r°), 

(a) Set£ m = (ti +t/)/2 and get r°(t m ) from r(t m ). 

(b) If r°(t m ) = r°(ii), set U = t m else set tf = t m . 

(c) Go on repeating steps (a) and (b) until tf — 
k = 1. 

2. Change n to n + 1 and go back to 1. 

Next we describe the method that we have adopted 
for restricting the trajectory of the state point of the 
system to a part of the configuration space. This is an 
extension of the interval bisection method described in 
the preceding paragraph. The trajectory is constructed 
part by part successively and in each part, whenever we 
detect a transition (using the interval bisection method) 
from the allowed part to forbidden part (say the part 
containing the basins of attraction of the four lowest- 
lying minima with V a < —43 in the case of the 13-atom 
Morse cluster), we re- initiate the MD at the point when 
the system last visited the allowed part using new ve- 
locities for the particles. We assign random velocities 
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drawn from a Maxwell distribution at temperature T [in 
the case of microcanonical simulation at constant £? tot , 
the temperature is obtained from the total kinetic en- 
ergy at the instant when the system last visited the al- 
lowed part i.e. T = 2(E tot - V)/(3JV - 6)]. After each 
such re-initiation we correct for the integrals of motion 
i.e. energy, momentum and angular momentum by suit- 
able rescaling and rotation of the velocities. 

Under this procedure the dynamics no longer remains 
truly Newtonian due to repeated velocity re-initiations. 
We have calculated the velocity autocorrelation function 
which decays very fast (within 2000 MD steps i.e. ~l-2 
vibrational periods) and heuristically we can argue that 
successive velocity re-initiations are more or less uncor- 
rected. Nonetheless, following reference [27] we have 
performed a few standard checks by computing quanti- 
ties that can be accessed through the above mentioned 
confined MD as well as through the regular unrestricted 
MD procedure and we found good agreement between 
the results obtained in these two different ways. For in- 
stance, the distributions of the total kinetic energy and 
the waiting time for a restricted part of the configura- 
tion space (e.g. the basin of attraction of a particular 
minimum) can be estimated [27] both by enforcing the 



trajectory to sample only the specified part or by picking 
out from a conventional long MD trajectory the intervals 
during which the system samples the specified part. 

We find that there are some spurious effects in the dy- 
namics due to cycling of trajectories near the saddles or 
the system spending more time near the border of the 
allowed and forbidden regions. Also, at low tempera- 
tures, when the system is constrained in the high energy 
parts such as Af r , it has a tendency to visit the lowest ly- 
ing basins more often. In the above mentioned procedure 
the system can visit the forbidden region within the small 
interval At and come back to the allowed part. We find 
these events to affect the properties related to short time 
relaxation such as the diffusion constant [Eq.(9)], waiting 
time t w [Eq.(ll)] and its distribution. While restricting 
the system in J\f r , as the system frequently escapes to 
the four lower lying basins, the barrier that appears in 
the temperature dependence of D and t w turns out to 
be related to the barriers connecting Af r to the forbidden 
low-energy part. However, since these visits are really 
short lived (< 100 - 200 MD steps), they do not affect 
the features of long-time relaxation, such as the decay of 
correlations characterized by Cg(t). 
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